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In this paper we use an approximation sequence defined by the 
Widder Laplace transform inversion formula to provide a practical 
method for inverting the Laplace transform. The approximation 
sequence converges uniformly and retains essential structural char- 
acteristics of the original function, e.g., nonnegativity, monotonicity, 
and convexity. Thus, we approximate a distribution function by 
distribution functions and use enhancement techniques to increase 
the speed of convergence and to capture the quality of exponential 
decay. Also, we present a practical computational method illustrated 
by examples. 

I. INTRODUCTION 

The purpose of this paper is to summarize the techniques presented 
in the paper "An Inversion Technique for the Laplace Transform" 1 
and to make available a useful reference of properties of the 
'approximation sequence,' and a new numerical method developed 
since the publication of Ref. 1. The inversion, or approximation se- 
quence, retains the essential structural characteristics of the original 
function, e.g., nonnegativity, monotonicity, and convexity. Thus, we 
approximate a distribution function by distribution functions. For 
application to queueing theory, this may be considered quite impor- 
tant. The basic inversion sequence, together with error estimates, is 
discussed in Section II; also, two enhancement procedures are given — 
namely, the construction of a sequence that is more rapidly convergent 
than the approximation sequence and which was not given in Ref. 1, 
and a method of accurate approximation to functions that decay 
exponentially. Section III discusses the new numerical method, and 
Section IV presents two examples of numerical inversion along with 
controls. Except for the new material whose derivations are given here, 
all proofs can be found in Ref. 1. 
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II. INVERSION SEQUENCE 

Consider a function f(t) for which the Laplace transform, f(s), 
defined by 

fta) = I e-«f(t)dt (1) 



-f-> 

Jo 



exists for s > c and f(t) = O(e c 0(£— > °°); then a sequence of functional, 
(L n )o, is defined by 



utn-tm-^*"** 



One has 



n+l" 



(2) 



lim f*(t)-f(t) (3) 

n— »oo 

uniformly in every finite closed interval throughout which f(t) is 
continuous. Equations (2) and (3) constitute a variant of Widder's 
theorem. 2 The sequence [f n (t)]o is called the "approximation se- 
quence" of f(t). 

The following lists some important properties of f n (t) valid for t > 0, 
n >: 0; a dot indicates differentiation. Assumptions on f(t) are to hold 
in (0, oo). 

(i) a < f(t) < b <=> a < /„(*) < 6, 

A(0+)=/(0+), 

/n(0+)=/(0+), 
/n(00)=/-(00), 

/„(oo)=/(oo), 

(ii) f(t) monotone => f n (t) monotone in the same sense. f(t) com- 
pletely monotone, absolutely monotone, convex, log-convex implies 
the same property, respectively, for f n (t). 

(Hi) f(t) > =» — — - monotone decreasing in n, f(t) convex «* 
n + 1 

f(t) < f n (t), f(t) > =* f n (t) monotone decreasing in n. 

(iv) If f(t)*g(t) = I fl-)g(u) — (Mellin convolution), then 

f n (t) is the approximation sequence of f(t) =* f n (t)*g(t) is the approx- 
imation sequence of f(t)*g(t). 
The pointwise error, €„(£; /), is defined by 

*n(t,f)=fn(t)-f(t). (4) 
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One has 



€»(£ f) ~ 2 TT G,(-n - 1, -n - l)t l f {l) (t), 
1=1 I- 



n 



(5) 



in which Gi(—n — 1, — n — 1) are Poisson-Charlier polynomials. 3 Set 
fit = Gi(—n — 1, —n — 1); then the following recursion is satisfied: 

fri— ^-r[A + A-!l ^o=l, 01 = 0. (6) 

n + l 

The terms of eq. (5) through order n~ 2 are 



2n+l 



+ 



(/i + D 2 
Some bounds for €„(£; /") are 



j*V-«*llTT«V-« 



/I 






for pointwise error, and 



| €„(*;/•) | <a„sup|*/(f) 

f>0 



uniform over all t. The a„ satisfies 

iin' 



OLn ~ 



00. 



(7) 



(8) 



(9) 



(10) 



It appears, therefore, that the convergence of f n (t) to f(t) is not 
rapid. In fact, from eq. (7) it is 0(l/n) and is not improved by assuming 
a higher degree of smoothness for f(t). However, one may decrease the 
error, e n {t, f), for a given n in at least two ways by: (i) modifying the 
sequence f n (t) to obtain a higher convergence rate, and (ii) modifying 
fo(t) to improve its approximation to f(t). The first way may destroy 
the desirable properties of the approximation sequence as previously 
listed because the transformation from f(t) to the members of the new 
sequence may no longer be positive; the second way can be imple- 
mented to preserve all the desirable properties with respect to a 
function different from the original f(t). 

The following procedure improves the convergence rate but does 
not preserve positivity. Define the functionals, o„, by 

o n f- o n (t) = (2 + ij Mt) - (1 + ij fn(t); (11) 

then, from eq. (7) 
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o n (t) - f(t) ~ - \ ) —r tY»(t) 

3 (n + l)(2n + 1) 

1 2n 2 + 9/i + 5 _,.,, , 

-8 ( n + l)'( 2n + l)^ y,(t) ' »— ■ (12 > 

and hence 

a n (t)-f(t)=0(l/n 2 ). (13) 

Often, fit) decreases exponentially fast for t — > oo; the inversion 
method may not capture this effect well for large t if the decrease is 
very rapid. This problem is largely overcome by the following tech- 
nique, which constitutes an example of the second method of improving 
the approximation. 

If f is) converges for s > -/?(/? > 0), and < a < /?, then f \s - a) 
exists and is the transform of git) = e at fit). If gnit) is the approxima- 
tion sequence for e at fit) and 

fnM) = e- at g n it), (14) 

then 

hmf n Jt)=fit). (15) 

71— *00 

The error e n , a it; f) = f n , a it) — fit) can be much smaller than e„(£ /), 
corresponding to a = 0, so that a real improvement in accuracy can 
result; of course, 

e n Jt;f) = e- at € n it;g). (16) 

The approximation f n , a it) imitates the rapid exponential decrease of 
fit) with increasing accuracy as a approaches /?. Since gnit) is an 
approximation sequence, it possesses all of the properties given earlier, 
many of which carry over to fit). 

III. NUMERICAL EVALUATION 

A type of generating function may be constructed for [/n(0]o. Let 

^v-iK 1 ! 1 ) (17) 



and define [a„(£)]o by 



then 



Giz,t)= I aM)z n ; (18) 

71=0 



fnit) = a n [-^j). (19) 
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This permits the use of numerical procedures for obtaining coeffi- 
cients in a power series. The following is such a method. The notation 

2irx 

e q (x) = e'^~ (20) 

will be used. Let < r < 1, q prime and q>n, and define the functional 
Sby 

Sf= ^±i J eq{ - nJ) f\l±l [i - re q (j)]}, (21) 

tqr y=1 [ t J 

then, in terms of f n (t), one has 

Sf= Ut) + £ /U (" + ^ 1 A r*. (22) 

For the purpose of computing f n (t),Sf in the form of eq. (21) will be 
used. 

We will now assume that f(t) is bounded so one may take 
\f(t) | < A. Clearly, from eq. (22) and property 1, 

\Sf-f n (t)\<Aj^ (23) 

and 

KmSf=fn(t). (24) 

r-*0+ 

Thus, Sf is an accurate approximation to f n (t) when r is small; 
however, if the round-off error of each term of eq. (21) is bounded by 
€, and the total round-off error by tj, then 

7,<e^-^r-\ (25) 

It follows that a choice of r may be made to ensure the round-off 
error, tj, is not too large, in accordance with eq. (25). The parameter q 
may now be chosen to render the truncation error given by eq. (23) 
comparably small. 

IV. EXAMPLES 

A Fortran program for the evaluation of Sf, as shown in Eq. (21), 
was written by Bharat Doshi. Double precision arithmetic was used 
with a round-off error of 10" 15 . As a rough estimate, it was assumed 
that the computations required to form each term of the summation 
in eq. (21) resulted in a round-off error of 10" 11 ; thus the choice e = 
10" 11 was made. The choices 
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r = 0.83, q = 127; n = 50, 
r = 0.91, q - 251; n = 100 



(26) 



lead to about five units error in the sixth decimal; the values of q given 
result in negligible truncation errors. These errors in computing fw and 
fioo were considered acceptably small. 
The function 



fit) = l/2e _l + l/2e- 



(27) 



is chosen for the first example. It represents a complementary distri- 
bution function whose values may be accurately computed as a control, 
and which shows the enhancement of accuracy obtainable by use of 
eqs. (11) and (14). Tables I and II illustrate these results. Since 

1 ' ' (28) 



hs)=l 



1 
+ - 



2s+l 2s + 3' 

the choice a = 1 is the appropriate value to use. 

It may be observed, from Table I, that the error is halved in going 
from /so to /ioo as expected from the 0(l/n) behavior of eq. (7), while 
the decrease of error from /ioo to am is better than 1/50 as seen from 
the 0(l/n 2 ) rate of eq. (13). Table II shows that the effect of a is 
greater for larger t since the exponential term that is being tracked 
becomes dominant. 

The second example concerns the complementary busy-period dis- 
tribution for an M/M/l. 4 The arrival rate is p and the service rate is 
/x = 1. For this case 

fit) = 1 - p" 1/2 \ e- <1+p) */i(2*v£) — , (29) 

Jo x 



and 



7(8) - ^" [P - 1 - 8 + V(p + 1 + S) 2 - 4p] , 

2ps 



(30) 







Table 1— 


■a 


= 




t 


f 


fm 




/ioo 


obq 


2 

4 
8 


0.06891 
0.00916 
0.00017 


0.07203 
0.01064 
0.00030 




0.07048 
0.00990 
0.00023 


0.06890 
0.00915 
0.00016 






Table II- 


-a 


= 1 




t 


f 


150.1 




/lOO.l 


OBO.l 


2 

4 
8 


0.06891 
0.00916 
0.00017 


0.06911 
0.00916 
0.00017 




0.06901 
0.00916 
0.00017 


0.06891 
0.00916 
0.00017 
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Table III— p = 


0.9, a = 




t 


f 


fa 


/ioo 


050 


3 
6 
9 


0.2919 
0.1950 
0.1514 


0.2942 
0.1967 
0.1528 


0.2931 
0.1959 
0.1521 


0.2919 
0.1950 
0.1514 




Table IV— p = 


0.5, a = 




t 


f 


fa 


/ioo 


"so 


3 
6 
9 


0.1803 
0.0764 
0.0399 


0.1834 
0.0785 
0.0415 


0.1819 
0.0775 
0.0407 


0.1803 
0.0764 
0.0399 




Table V— p = 


0.1, a = 




t 


/ 


fa 


/ioo 


CT50 


3 
6 
9 


0.0732 
0.0088 
0.0013 


0.0774 
0.0103 
0.0018 


0.0753 
0.0096 
0.0016 


0.0732 
0.0088 
0.0013 




Table VI- 


-p « 0.1, a* 0.4675 


* 


f 


fa. 


/lOO,a 


O50,a 


3 

6 
9 


0.0732 
0.0088 
0.0013 


0.0742 
0.0090 
0.0014 


0.0737 
0.0089 
0.0014 


0.0732 
0.0088 
0.0013 



in which Ji(2*vp) is a modified Bessel function. The evaluation of the 
control was accomplished by using a quadrature procedure. The values 
p = 0.9, 0.5, and 0.1 were used. The best choice of a is, from Section II, 
(1 — Vp) 2 ; however, when a is small, one could simply set a = 0, since 
there is little improvement in using the best a. For the case p = 0.1, 
one has a = 0.4675; accordingly, a comparison was made with a = 0. 
Tables III through VI sample the results obtained. 

V. CONCLUSION 

It may be observed that the best a for p = 0.1 noticeably improved 
the approximation. The performance of the o„ functionals in all cases 
created a marked improvement in accuracy; however, one should 
remember that the desirable properties of the L n functionals are not 
possessed by the a„. Nonetheless, there is a property of the a n , impor- 
tant in numerical work, which follows from eq. (11). Namely, if 8 is the 
round-off error in the computation of f n (t) and fon(t), then (3 + 2/n)8, 
n > 1 bounds the induced round-off error in ct„ (i.e., round-off errors 
are not significantly magnified). 
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An interesting use of the inversion technique follows. A function f(t) 
is given in the form 



fit) = g(t, x)dx, (31) 

Jo 

with g(t, x) known; it is required to evaluate f(t) over a wide range of 
values of t. Quadrature methods are accurate, for t small, or can be 
designed, for t large, but do not apply equally well over all values of t, 
including the important transition region from small to large. If the 
transform f(s) can be obtained, then the inversion method described 
here can be used to obtain sufficiently accurate values of f(t) over the 
entire range of t, since eq. (9) shows that the errors, € n (t; f), are 
uniformly bounded. 
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